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Summary: 


The  simply  coupled  scale  analysis  method  (SCSAM)  has  been  further  extended  to 
perform  fracture  analysis  of  heterogeneous  materials  and  multilayered  composite 
materials  in  particular.  This  work  is  based  on  Breitzman,  et  al.  [1]  and  Lipton  [2,3] 
SCSAM  work  allowing  prediction  of  the  local  stress  and  strain  fields  inside  undamaged 
pre-stressed  composites.  SCSAM  allowed  for  rigorous  micro-level  stress  bound 
evaluation  in  the  presence  of  stress  gradients,  including  regions  where  periodicity  is  not 
maintained  in  all  directions.  This  methodology  allowed  evaluating  the  strength  of 
composite  scarf  repair  and  optimizing  the  repair  patch  stacking  sequence  to  design  repair 
satisfying  low  observability  requirements.  Although  basic  results  on  repair  analysis  were 
reported  in  the  previous  FY04-06  period,  the  results  of  the  optimization  study  were 
obtained  later  and  are  included  in  the  present  report.  Critical  development  of  the  present 
project  is  an  extension  of  SCSAM  to  simulate  the  failure  process,  whereas  previously  wc 
evaluated  strength  based  on  stress  levels  in  undamaged  composites.  Progressive  multi- 
eraek  fracture  simulation  on  each  scale  is  performed  by  using  the  mesh-independent 
crack  propagation  technique  and  cohesive  formulation  for  crack  extension.  A 
representative  volume  element  (RVE)  approach  employed  for  homogenized  stiffness  and 
stress  concentration  factor  computation  was  extended  to  compute  homogenized  cohesive 
crack  opening  relationship  on  the  higher  hierarchical  scale  based  on  the  cohesive 
relationships  of  the  constituents.  In  the  spirit  of  SCSAM  and  within  the  limitation  of 
sufficiently  small  microstructure  scale  simulation,  the  coupling  of  scales  is  occurring 
through  hierarchical  computation  of  homogenized  properties,  which  now  include  fracture 
properties,  whereas  the  simulation  of  the  failure  process  on  eaeh  scale  is  eondueted 
independently.  Such  an  approach  is  highly  practical  and  applicable  to  prediction  of 
complex  phenomena  in  composite  and  hybrid  materials.  An  example  of  application  of 
this  methodology  is  simulation  of  thermo-oxidative  cracking  in  high-temperature 
composites,  which  we  eondueted  recently  in  support  of  the  Hybrids  for  Extreme 
Environments  Program  in  AFRL/RXBC;  Marilyn  Unroe  is  the  program  manager,  c-mail; 
Marilyn. Unrocrdjwpafb.af  mil  and  Dr.  Vernon  Beehel  technical  lead,  c-mail: 
Vernon. BcchcKojwpafb.af  mil.  We  have  also  performed  simulations  of  complex  failure 
patterns  in  multilayered  composites  under  more  traditional  loading  in  ambient  conditions. 

Further  research  is  required  to  address  the  anisotropy  of  the  fracture  toughness  in 
materials  with  strongly  anisotropic  homogenized  properties  and  related  questions  of  the 
influence  of  the  RVE  on  the  results  of  the  homogenized  cohesive  relationship  prediction. 
A  critical  underlying  technology  for  application  and  development  of  progressive  fracture 
modeling,  SCSAM  is  a  technique  for  modeling  complex  multiple  site  initiation  and  eraek 
propagation  patterns  in  three-dimensional  formulation,  which  also  requires  further 
development.  It  is  also  envisioned  that  the  future  development  path  of  the  SCSAM 
framework  should  include  nonclassical  crack  models  at  a  small  length  scale. 
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1.  Introduction  and  State  of  the  Art 

To  date,  the  primary  approach  used  for  multi-scale  failure  prediction  in  composite 
material  is  the  so-ealled  continuum  damage  mechanics  (CDM)  method  [4-6].  According 
to  this  method  the  effect  of  damage  at  each  length  scale  is  manifested  by  local  stiffness 
degradation  at  the  higher  (coarser)  length  scale.  The  attractiveness  of  this  technique  is  in 
the  intuitiveness  of  its  concept  and  the  ease  of  implementation.  One  reason  for  the 
popularity  of  CDM  is  the  potential  ability  to  model  an  arbitrary  path  of  damage 
propagation  on  the  macro-scale  using  numerical  schemes  based  upon  a  standard  finite 
element  (FE)  formulation.  The  drawback  of  these  models  is  their  inability  to  accurately 
describe  the  local  effects  of  interaction  of  various  damage  modes  and  local  effects  of 
stress  redistribution  in  the  damaged  area  as  discussed  by  larve,  et  al.  [7].  Instead,  a 
computational  framework  based  on  discrete  modeling  of  arbitrary  multiple  interacting 
cracking  is  chosen  for  development  and  implementation  of  SCSAM. 

Several  methodologies  have  been  proposed  for  modeling  the  kinematics  of 
arbitrary  cracking.  Evolution  of  a  crack  front  can  be  captured  by  traditional  FE  modeling 
combined  with  adaptive  remeshing  techniques  [8].  Successful  in  predicting  complex 
crack  evolution  in  metallic  structures,  its  application  to  predicting  crack  scenarios  in 
laminated  composites,  where  cracks  form  in  different  plies  in  adjacent  locations,  will 
require  remeshing  of  significant  volumes  under  multiple  mesh  compatibility  constraints. 
An  alternative  approach  to  modeling  crack-induccd  displacement  discontinuities  involves 
mesh  independent  crack  modeling  techniques.  Early  works  devoted  to  mesh  independent 
modeling  of  matrix  cracking  in  composite  laminates  include  [9,  10].  Over  the  last 
decade,  following  the  pioneering  work  of  Moes,  et  al.  [11]  in  which  the  concept  of  the 
extended  Finite  Element  Method  (x-FEM)  was  proposed,  a  significant  effort  to  further 
develop  these  ideas  has  been  put  forth.  Although  most  of  the  research  was  devoted  to 
arbitrary  crack  propagation  in  isotropic  materials,  recent  application  to  composite 
materials  includes  delamination  modeling  and  textile  composite  architecture 
representation  [12,  13]. 

In  the  present  development  we  target  complex  materials  where  several 
hierarchical  scales  exist.  Two  types  of  cracking  need  to  be  generally  recognized  in  such 
systems:  the  interface  cracking  and  arbitrary  three-dimensional  intra-phase  cracking,  as 
shown  in  Figure  1.  This  classification  is  applicable  on  each  hierarchical  scale,  i.e.  Figure 
1  shows  a  schematic  view  of  a  multilayered  composite  skin.  Figure  l.b  shows  both 
interface  and  intra-phase  (ply)  cracking  on  the  ply  scale  level.  Interface  and  intra-phase 
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Figure  1.  Hierarchical  composite  structure 


matric  cracking  on  the  fiber  matrix  level  is  shown  in  Figure  1  .c.  In  this  case  the  interface 
cracking  is  between  the  matrix  and  the  individual  fibers,  and  the  intra-phase  cracking  is 
inside  the  matrix.  In  order  to  facilitate  the  SCSAM,  a  computational  framework  capable 
of  originating  and  evolving  complex  interface  and  intra-phase  cracking  networks  without 
any  prior  knowledge  or  assumptions  regarding  the  locations  of  damage  origination  is 
required. 

Our  analysis  is  built  around  a  regularized  x-FEM  formulation  (larve  [14]),  which 
is  extended  to  multiple  crack  initiation  and  evolution  modeling.  In  a  regularized 
formulation,  the  step  function  describing  the  crack  surface  is  replaced  by  a  continuous 
function  and  consequently  allows  one  to  maintain  a  Gauss  integration  schema  for  element 
stiffness  matrix  computation  without  regard  to  crack  orientation.  In  this  case  the 
connection  between  two  arbitrarily  cracked  surfaces  can  be  easily  established  by 
computing  integrals  of  the  products  of  the  shape  functions  on  the  two  surfaces. 

We  begin  by  outlining  the  regularized  mesh  independent  crack  (MIC)  modeling 
framework,  followed  by  the  nonlinear  homogenization  of  the  RVE  containing  evolving 
cracks.  Next  we  will  discuss  computation  of  the  homogenized  effective  stiffness  and 
strength  and  fracture  toughness  properties  of  high-temperature  unidirectional  G3- 
500/PMR15  composites  subjected  to  thermal  oxidation.  Finally,  we  illustrate  complex 
fracture  phenomena  in  a  laminated  composite  such  as  is  shown  in  Figure  1 . 


2.  Mesh  Independent  Intra-Phase  and  Interface  Crack  Modeling 


As  mentioned  before  we  seek  a  framework  capable  of  accommodating  the 
initiation  and  interactive  evolution  of  both  intra-phase  and  interface  cracking.  Without 
restricting  the  generality,  we  shall  consider  a  laminated  media,  where  the  interface  cracks 
represent  delamination  and  the  intra-phase  cracks  are  transverse  cracks.  An  application  to 
arbitrary  surface  geometry,  where  the  interface  cracks  will  represent  fiber  surfaces,  will 
be  considered  in  the  next  section. 

a.  Mesh  Independent  Crack  Modeling 

Consider  a  partition  of  unity  set  of  continuous  basis  functions  X,{xJ  and  a 
displacement  approximation  on  the  domain  of  interest  F 
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(1) 


u(x)  =  ^X,(x)U, 

/en 


Next  consider  a  crack  appearing  in  this  volume,  as  shown  in  Figure  2,  with  a  surface  Fa 
defined,  by  means  of  the  signed  distance  function  defined  as: 


fa  (x)  =  •s/g«(n(x)(x  -  x))  min  X  -  x 


(2) 


where  n(x)  is  the  normal  to  the  crack  surface  f  at  the  point  x.  The  traditional  x-FEM 
strong  discontinuity  formulation  is  based  on  the  element  enrichment  with  displacement 
modes  discontinuous  over  the  crack  surface  by  introducing  shape  functions  multiplied  by 
HifJji)),  where  H(x)  is  the  Heaviside  step  function.  Consider  a  partition  of  unity  set  of 
continuous  basis  functions  Xi(\)  and  a  displacement  approximation  on  the  domain  of 
interest  V. 


Figure  2.  Signed  distance  function  for  arbitrary  matrix  crack  definition  and  schematics 
of  cross  integral  computation  for  adjacent  ply  interfaces  in  the  regularized 
formulation. 
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In  the  regularized  formulation  [14]  the  Heaviside  step  function  is  replaced  with  a 
continuous  function  H{\) 


H{x)  =  ^X,{x)h, 

/en 


(3) 


where  Xi(x)  are  the  same  shape  functions  as  in  Eqn.  (  1).  The  coefficients  hi  are  calculated 
as  follows 


l/  I  fy  /«(x)A^i(x)dx\ 

Jv,  l/«(x)|A'i(x)dxy  (4) 


This  definition  involves  only  continuous  functions  and  can  be  calculated  by  using 
standard  Gauss  quadratures.  The  only  coefficients  h,-,  which  are  not  equal  to  0  or  1, 
correspond  to  those  shape  functions  that  are  divided  by  the  crack  surface.  Let  us  denote 
the  set  of  such  index  values  (for  which  hi  is  not  equal  to  0  or  1)  by  Qq.  Then  the  enriched 
approximation  for  the  domain  V  and  arbitrary  crack  is  defined  in  the  following  form 


and 


u  =  //u"’+(l-//)u''’  +  u‘'’ 


/€n/Q^ 


(5) 

(6) 


(7) 


We  have  also  omitted  the  spatial  argument  for  conciseness.  The  displacement 
approximation  given  in  Eqn.  (9)  contains  the  enrichment  in  the  crack  region  via 
displacement  fields  and  as  well  as  the  unchanged  displacement  field,  away 
from  the  crack.  Equations  (  5)-(  7)  define  the  enriched  displacement  approximation  by 
replacing  each  original  shape  function  JG  influenced  by  the  crack,  /E  Qa,  with  two  shape 

functions,  H  Xi  and  (1-//)A/.  This  approximation  was  applied  in  [14]  in  conjunction 
with  a  higher  order  displacement  approximation  (p-elcments)  as  well  as  B-splinc 
approximation  of  displacements,  where  the  coefficients  Ui  do  not  correspond  to  nodal 
displacements.  The  convenience  of  such  representation  of  the  enrichment  lies  in  the 
simplicity  of  bookkeeping  the  connectivity,  where  the  two  copies  of  the  shape  function 

do  not  interact  and  are  connected  with  only  alike  H  or  (1-^)  multiple  copies  of  the 
other  enriched  shape  function  if  their  supports  overlap. 
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The  strain  energy  in  the  volume  V  with  displacement  approximation  given  in  Eqn. 
(  5)  is  written  as 


^  —  j  -  e)  +  (l  —  —  e)  ^ 

+  —  e)  +  —  e)^C(£^^^  —  e) 

+  (l  —  —  e)^C(£^^^  —  e)}dV 


where  superscript  T  means  transpose  operation,  C  is  the  elastic  orthotropie  stiffness 
tensor,  and  e  is  the  nonmechanical  strain  tensor  introduced  to  account  for  thermal 
processing  stresses,  so  that  e=(T-To)a.  Here  T  and  To  are  cure  and  ambient  temperature, 
and  a  is  the  thermal  expansion  tensor.  The  components  of  the  stiffness  tensor  in  the 
arbitrary  coordinate  system  can  expressed  according  to  [15].  The  strain  tensors  are 
computed  from  displacement  fields  in  Eqn.  (  5),  so  that 

fW  =  1  (VuW  +  ^ 


Expression  (  9)  does  not  contain  the  interaction  energy  in  the  gradient  zone,  which 
will  be  considered  next.  Typically,  in  a  regularized  formulation,  the  MIC  propagation  is 
governed  by  the  constitutive  properties  in  the  gradient  zone,  in  which  |V//|  = 

yJVH  •  VH  >  0.  In  the  present  formulation,  the  cohesive  constitutive  relationship  [16] 
developed  for  interface  fracture  modeling  is  inserted  into  the  gradient  zone  of  the 
regularized  formulation  directly.  It  can  be  accomplished  by  using  Dirac’s  delta  function 
of  signed  distance  function  of  the  crack  surface  to  express  the  cohesive  energy  of  the 
crack  surface  through  the  volume  integral  in  the  gradient  zone.  Indeed  the  surface  with 
signed  distance  function  enclosed  in  an  arbitrary  volume  v  can  be  expressed  as  a 
volume  integral 


s.  =  lsM)dv  (10) 

V 

where  ^£,(7^)  is  the  Dirac’s  delta  function  of  signed  distance  function.  One  can  also 
establish  that  for  an  arbitrary  continuous  function  g(x)  defined  in  volume  V,  a 
relationship  between  the  surface  integral  over  the  crack  surface  Fa  (FacV)  and  a  volume 
integral  exists  so  that 


||g(x)r/5  =  |g(x)^o(/«V^  (11) 
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This  relationship  can  be  readily  established  by  applying  Eqn.  (  10)  in  small  adjoining 
volumes  encompassing  surface  Fa  to  develop  the  integral  sums  representing  the  left-  and 
right-hand  side  of  Eqn.  (  11).  In  the  case  of  the  regularized  formulation  we  compute  the 
approximate  value  of  the  right-hand  side  by  replacing  the  Dirac’s  function  ^^(/a)  the 

gradient  of  the  approximate  step  funetion|V//|.  The  eontinuous  function  g(x)  defined  over 
the  volume  is  now  replaeed  with  the  point  wise  eohesive  energy  of  craek  opening.  In  the 
traditional  interfaee  eohesive  formulation,  the  cohesive  energy  is  a  funetion  of  the  eraek 
opening  displacement  and  is  also  dependent  on  the  ratio  of  the  opening  mode  I 
(perpendicular  to  erack  surface)  and  the  shearing  mode  II  tangential  to  the  erack  surface. 
To  separate  these  modes  one  needs  to  know  the  direetion  of  the  normal  to  the  eraek 
surfaee  in  all  points.  In  the  regularized  formulation,  the  displaeemcnt  gap  and  the  normal 
vector  arc  defined  in  all  points  of  the  gradient  zone  as 


and  n  =  V/« 


(12) 


where  and  M^^^are  defined  by  using  the  displaeement  fields  in  enriehed 
approximation  (  5).  The  eohesive  energy  of  the  eraek  opening  g(x)  will  be  eonsidered 
homogeneous  and  therefore  dependent  upon  the  spatial  eoordinate  x  only  as  a  funetion  of 
the  displaeement  gap  and  the  normal  veetor  to  the  crack  surfaee,  so  that  g(x)=g(AM,  n). 
Thus  the  fracture  energy  of  the  MIC  propagation  is  expressed  as 

M  =  ||VW|g(AM,n)dV.  (13) 

V 

b.  Cohesive  Energy  and  Interface  Crack  Modeling 

The  shape  of  the  interface  cohesive  energy  funetion  in  composite  laminates  for 
predicting  the  delamination  failure  mode  has  being  studied  extensively  ineluding  [16-18]. 
It  ean  be  written  down  in  the  invariant  form  a  function  of  the  absolute  value  of  the 
displacement  gap  ^=|A«|  and  a  mode  mixity  parameter,  such  as 

^  ^ - - ’ 

introduced  in  [1 6],  where  Au„  =  (Am  •  n)  is  the  normal  component  of  the  displacement 
gap,  which  is  equal  to  0  for  mode  I  propagation  and  I  for  mode  II  propagation.  The  shape 
of  the  relationship  g(^,B)  ean  vary;  however,  to  assure  the  correct  crack  propagation 
eharaeteristies,  it  must  satisfy  the  following  condition 


g(c«,B)=Gc(B), 


(  14) 


Gc(B)  =  G„  +  (G„  -  G„c)B\ 


lO 


where  G/c  and  Guc  are  experimentally  measured  fraeture  toughness  values  and  ti=2.25. 
The  functional  shape  of  the  fraeture  energy  as  a  function  of  the  displacement  gap  is 
defined  by  the  relationship  between  the  cohesive  tractions  and  the  displacement  gap, 
which  in  [16]  is  assumed: 


r  =  (l-d)KAu  +  dK{Au-n),  (15) 

where  K  is  high  initial  stiffness  and  d  is  the  damage  parameter.  The  first  term  in  Eqn. 
(15)  represents  the  eraek  cohesive  force,  and  the  second  term  prevents  interpenetration. 

The  brackets  {x)  =  ^  (x  +  |x:|)  represent  the  MeAuley  operator.  A  bilinear  relationship  is 

assumed  for  T(k)  such  that  d=0  if  Ao  and  d=l  if  Af .  The  cohesive  energy  g(^,B)  is 
the  area  under  the  T(k)  curve,  and  thereby  according  to  Eqn.(  14),  the  final  value  of  the 
displacement  gap  is  determined  by  the  its  initial  value  and  the  fraeture  toughness  as 

Ar2Gc/(KAo) 

The  initial  value  of  the  displacement  gap,  beyond  which  the  interface  failure  begins,  is 
defined  as 


Ao=to/K, 


where  Vg  is  the  initial  strength,  also  depending  on  the  on  the  mode  mixity  parameter  B  as 
follows 


(joY  =  Y^  +  (S^  -  Y^)B^ 

where  Y  and  S  are  the  transverse  and  shear  strengths,  respectively.  All  parameters 
entering  the  analysis,  such  as  the  fraeture  toughness  and  strength  values,  represent 
material  properties  and  are  measured  by  using  standard  test  methods.  The  reader  is 
referred  to  reference  [16]  for  full  details,  and  the  brief  description  above  is  given  for 
completeness  of  the  present  formulation. 

The  cohesive  energy  associated  with  the  delamination  of  the  plies  is  computed  by 
integrating  the  cohesive  energy  over  the  interface  between  plies  n  and  n+1,  separated  by 
the  horizontal  surface  z=Zn  is 


JJg(X.B)ds  (16) 

z=Zn 

In  this  ease  the  normal  veetor  to  the  eraek  surfaee  is  (0,0,1),  and  the  displacement  gap 
between  the  interfaces  is  computed  by  using  enriched  displacement  approximations  (  5) 
in  all  plies. 


c.  Variational  Formulation  and  Problem  Solution 


The  work  of  the  external  traetion  P  applied  at  the  portion  Sp  of  the  boundary  dV 
is  expressed  as 


/4=jPuc/S  (17) 

Sp 

To  derive  the  equilibrium  equations  in  terms  of  the  displaeement  approximation 
eoefficients  for  a  body  whieh  eontains  multiple  matrix  and  delamination  eraeks,  the  first 
variation  of  the  potential  energy  is  required  to  vanish.  Combining  Eqs. 

(  8),  (  13),  (  16)  and  (  17)  for  the  volume  strain  energy  and  the  work  of  the  surfaee 
tractions,  and  summing  over  all  plies  and  interfaces,  gives. 


<y(E^[=l(Wn  +  Mn  )  +  “  ^)=0 


(18) 


The  lower  indexes  in  strain  energy  Wn  and  the  MIC  cohesive  energy  Mn  designate  that 
they  arc  computed  for  the  n-th  ply;  the  delamination  free  energy  On  is  between  the  n  and 
n+l-th  ply.  Performing  the  variation,  the  following  system  of  equations  is  obtained 


fW+M-HD;f/=/»+Af  (  19) 

The  vector  of  unknowns  is  arranged  by  ply  in  the  order  U^=(Ul'l,Ul^l,ul^l,...)^.  The 
matrix  W  and  right-hand-side  vector  N  are  obtained  from  a  variation  of  W„  and  are  the 
elastic  stiffness  matrix  and  the  nonmechanical  load  vector,  where  the  shape  functions 

and/or  their  derivatives  in  each  integration  point  are  multiplied  by  H  or  by  (1-//),  or 
unmodified. 

In  the  case  when  MlCs  are  present  in  adjacent  plies,  i.e.  n  and  n+1,  the 
displacements  Un  and  Un+i  arc  enriched  according  to  Eqn.  (  5).  When  regularized 
enrichment  is  used,  the  integration  of  all  shape  functions  is  conducted  on  the  original 
Gauss  points,  and  thereby  is  straightforward  for  arbitrary  enrichment  of  displacements  in 
the  adjacent  plies.  Two  sparse  matrices  F{{  and  are  formed  at  the  interface  between 
plies  n  and  n+1.  For  each  integration  point  (always  common  for  plies  n  and  n+1)  on  the 
common  interface  (row  index),  these  matrices  contain  the  values  of  all  nonzero  shape 
functions  in  this  point  (column  index)  for  the  ply  n  and  ply  n+1,  respectively.  These 
matrices  have  large  dimensions  but  arc  very  sparse.  Matrix  O  can  then  be  represented  in 
the  following  form 

w-i 

^  (f -  Fnf  ^(^■r '  -  f H)  (  20) 

n=l 
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The  weight  matrix  Z)  is  a  diagonal  matrix  eontaining  the  Gauss  weights  and  surfaec 
Jacobian  at  the  given  integration  point  multiplied  by  the  cohesive  stiffness  (l-d)K.  Note 
that  the  matrices  F„  can  also  be  used  to  compute  the  displacement  gap  Au  in  each 
integration  point  by  simple  multiplication  —  F^)U.  By  knowing  the  displacement 

gap,  one  can  update  the  damage  variable  and  the  coefficients  of  decohcsive  stiffness  in 
each  integration  point. 

The  matrix  M  has  a  structure  similar  to  O.  For  opening  of  the  MICs,  the  matrices 
of  shape  function  values  in  integration  points  are  formed  so  that  F\  and  f  ^  contain  the 
functions  corresponding  to  enrichments  and  in  Eqn.  (  5),  respectively,  and  the 
components  of  the  weight  matrix,  which  wc  denote  D  to  distinguish  from  the 
delamination  case  D,  are  the  products  of  Gauss  weight,  Jacobian,  |V/j|  and  the 
decohesive  stiffness  (l-d)K  in  the  given  integration  point.  Thus  the  contribution  of 
variation  of  the  MIC  opening  energy  can  be  written  as 

N 

M  =  Y^{Fl-FiyD{Fl-Fl).  (21) 

n=l 

The  loading  vector  P  comprises  of  applied  external  stress  as  well  as  displacement 
loading.  Let  the  vector  of  displacement  approximation  coeffieients  at  the  edges  x=0,  L  of 
the  thermal  expansion  problem  with  the  boundary  conditions  be  Qo.  Then  the  load  vector 
at  the  load  step  k  is 

P=Qo+QrY.U^i 

where  0/  is  a  vector  of  unit  vector  displacement  at  x=0,L  and  Aj-  are  load  increments  at 
the  /-th  step. 

The  system  of  equations  (21)  contains  highly  nonlinear  components  M  and  O 
and  is  solved  by  using  the  Newton-Raphson  (NR)  method  at  each  step. 

3.  Intra-Phase  (Ply)  Mic  and  Interface  Cracking  (Delamination)  Interaction 

Verification  Study 

Numerical  results  devoted  to  verification  and  application  of  the  proposed 
methodology  for  failure  prediction  in  laminated  composites  will  be  presented  below.  The 
verification  will  consist  of  modeling  failure  in  specimens  with  predefined  matrix  cracking 
patterns  and  comparison  with  experimental  data  and  FE  analysis.  The  purpose  of  the 
following  example  is  to  test  the  methodology  of  interface  and  intra-phase  cracking, 
represented  by  MIC.  In  this  case  we  are  less  eoneemed  with  hierarchical  reproduction  of 
the  orthotropic  unidirectional  composite  ply  properties  entering  this  analysis.  However,  it 
is  pertinent  to  discuss  the  set  of  material  properties  required  for  a  ply  level  fracture 
simulation  in  a  composite  laminate.  These  properties  pertain  to  a  unidirectional  ply  and 
are  listed  for  the  material  system  T300/914C  in  Table  1.  We  have  divided  these 
properties  in  two  groups.  The  first  group  represents  stiffness  and  thermal  expansion 
characteristics,  which  can  be  readily  computed  by  using  methodology  given  by 
Breitzman,  et  al.  [1]  from  the  fiber  matrix  level  constituent.  It  is  the  second  group  of 
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properties,  homogenization  of  which  we  will  be  discussing  in  the  next  section.  At  present 
we  will  use  the  values  given  in  Table  1  obtained  experimentally. 


Table  1 .  Material  properties  used  for  validation  and  strength  prediction  examples 


Comments 

Property 

T300/914C 

IM7/8552 

Ref  [20] 

Ref  [32] 

En(GPa) 

139.9 

161 

Group  I ; 

E22,  E33  (GPa) 

10.1 

11.38 

Stiffness 

G23(GPa) 

3.7 

3.98 

and 

G,2,i3  (GPa) 

4.6 

5.17 

Thermal 

Expansion 

Properties 

V23 

0.436 

0.436 

Vl2,Vl3 

0.3 

0.32 

an  (  l/'C) 

0.4x  10'^ 

0 

Ot22,  Ct33  (  1/^C) 

2.5x  10'^ 

3x  10'^ 

Processing 

T-To  (®C) 

-150 

-150 

Group  II: 

Yt  (MPa) 

80 

60 

Strength 

Yc  (MPa) 

300 

260 

and 

S  (MPa) 

100 

90 

Fracture 

Toughness 

G,c(J/M') 

120 

200 

G„c(J/M^) 

500 

1000 

a.  Transverse  Crack  Tension  Test  (TCT) 

The  TCT  specimen  described  in  reference  [20]  was  examined  to  evaluate  the 
accuracy  of  delamination  propagation  prediction,  which  emanates  from  a  MIC.  This 
specimen  consists  of  three  unidirectional  plies  with  thickness  t,  2t  and  t,  respectively,  and 
is  subjected  to  axial  tension.  The  middle  ply  is  cut  at  the  center  through  the  width.  When 
the  tensile  loading  is  applied  at  the  x=0,L  edges  in  the  x-direction,  the  applied  stress  at 
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first  increases  linearly  until  the  crack  in  the  middle  ply  will  initiate  delaminations 
between  the  middle  and  the  outer  plies.  These  delaminations  will  then  propagate  in  a 
stable  manner  keeping  the  applied  stress  constant.  Finally  after  the  delaminations  reach 
the  grips,  the  load  will  start  increasing  again  but  with  a  different  slope.  This  configuration 
allows  one  to  examine  the  effect  of  delamination  origination  from  a  single  matrix  crack 
and  was  chosen  for  model  validation  purposes.  The  problem  was  solved  both  by  using  the 
standard  FE  model,  shown  in  Figure  3a,  and  the  MIC  model,  shown  in  Figure  3b.  In  both 
eases,  cohesive  interfaces  were  implemented  on  the  surfaces  between  plies.  Only  half  the 
laminate  was  modeled,  i.e.  two  plies  with  equal  thickness  t  and  zero  z-displaeement 
condition  on  the  surface  z=0.  In  the  ease  of  the  FE  model,  the  crack  in  the  middle  ply  is 
aligned  with  a  mesh  line  and  simply  modeled  by  using  double  nodes.  No  mesh  line 
aligned  with  the  initial  crack  exists  in  the  case  of  curved  nonuniform  mesh  in  Figure  3b, 
and  the  MIC  model  is  used  to  insert  the  eraek  in  the  middle  ply  only.  The  sign  distance 
function  is  /a=x-L/2.  The  total  number  of  axial  intervals  in  the  two  models  coincides; 
however,  due  to  mesh  nonuniformity,  the  local  density  of  the  MIC  mesh  varies 
significantly.  Applied  displacement  versus  applied  load,  predicted  by  the  two  models,  is 
shown  in  Figure  4  and  are  nearly  identical. 


Y 


(•) 


(b) 


Figure  3.  Transverse  eraek  tension  (TCT)  specimen.  Conventional  model  (a)  and  MIC 
evaluation  model  with  skewed  mesh. 
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Figure  4.  Load  displacement  curves  of  the  TCT  specimen  testing  predicted  by  using 
conventional  FE  model  and  MIC. 

4,  Effective  Fracture  Toughness  Computation 

Fracture  toughness  of  composite  materials  is  a  strongly  anisotropic  characteristic. 
References  [21,  22]  capture  the  development  of  fracture  toughness  measurement 
standards  and  the  physics  of  mixed-mode  delamination  in  unidirectional  composites. 
Very  little  work  has  been  performed  on  experimental  investigation  of  the  relationship 
between  fracture  toughness  of  the  matrix  and  fiber  matrix  interfaces  and  that  of  the 
effective  fracture  toughness  of  the  composite  [23,24].  Below  we  attempt  to  predict  these 
characteristics.  We  first  describe  the  simplest  RVE,  which  is  capable  of  representing  the 
stiffness  portion  of  the  composite  properties,  i.e.  the  Cubic  Unit  Cell,  and  on  its  example, 
the  fracture  toughness  homogenization  methodology.  As  will  be  shown,  the  minimum 
requirements  to  an  RVE  for  fracture  toughness  characterization  arc  not  satisfied  for  this 
configuration,  and  next  we  consider  the  Hexagonal  Unit  Cell.  We  end  with  the 
application  of  the  developed  methodology  to  Thermal  Oxidation,  and  Future  Directions. 
Only  Mode-I  fracture  toughness  of  a  unidirectional  composite  has  been  addressed  at 
present. 

The  Mode  I  fracture  toughness  (G[c)  of  a  composite  is  experimentally  determined 
from  the  double  cantilever  beam  (DCB)  test  schematically  shown  in  Figure  5.  At  a 
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critical  tensile  load  applied  to  the  ligaments,  a  pre-existing  crack  is  starting  to  advance. 
This  load  is  recorded  and  used  to  evaluate  the  fracture  toughness  at  the  crack  tip  at  the 
propagation  onset.  The  fracture  toughness  determined  from  this  experiment  is  defined  as 
the  Mode  1  fracture  toughness  of  the  composite. 


Figure  5.  DCB  test  schematic  for  Mode  1  fracture  toughness. 


The  purpose  of  this  work  is  to  understand  how  to  predict  Mode  1  fracture 
toughness  at  the  ply  length  scale  from  simulations  at  the  fiber/matrix  microstructural 
length  scale. 

a.  Cubic  Unit  Cell 

Perhaps  the  simplest  unit  cell  used  to  represent  a  fiber-rcinforced  composite  is  the 
cubic  unit  cell.  A  two-fiber  model  (see  Figure  6)  is  used  to  predict  ply-lcvel  Mode  1 
fracture  toughness  for  the  composite.  The  dimensions  of  the  two-fiber  model  are 
LxLx2L  (jcx^xz),  with  the  origin  of  the  coordinate  system  as  shown  in  Figure  6.  For 
illustrative  purposes,  red  arrows  have  been  added  to  show  the  location  of  the  cohesive 
zone  interface.  The  cohesive  zone  is  the  fracture  surface  for  this  model.  The  cohesive 
zone  interface  model  is  defined  across  an  interface  as  a  combination  of  normal 
displacement,  total  displacement,  surface  tractions,  and  fracture  toughness  as  previously 
described.  For  all  simulations  in  this  section,  a  60%  fiber  volume  fraction  is  assumed 
unless  otherwise  noted.  Input  stiffness  and  thermal  expansion  properties  for  the  fiber  and 
matrix  phases  are  summarized  in  Table  2  and  are  from  reference  [24].  We  also  show  the 
Group  I  stiffness  and  thermal  expansion  properties  computed  according  to  reference  [1]. 
In  addition  to  stiffness  and  thermal  expansion  properties,  the  fracture  simulations  require 
four  additional  input  properties.  These  properties  are  the  strength  and  fracture  toughness 
of  the  matrix  material  and  of  the  fiber/matrix  interface.  The  matrix  strength  and  fracture 
toughness  are  77  MPa  and  177  J/m^,  respectively.  The  interface  of  this  material  system  is 
chemically  optimized  and  is  thus  quite  strong.  Without  specific  data  on  the  correct 
intcrfacial  values,  we  set  the  interfacial  strength  and  fracture  toughness  equal  to  those  of 
the  matrix  material.  Wc  do  not  allow  the  fibers  to  fracture.  As  a  reference,  the 
experimentally-determined  composite  strength  and  Mode  1  fracture  toughness  values  to 
which  we  will  be  comparing  are  66  MPa  and  225  J/m^,  respectively. 
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Figure  6.  Cubic  unit  cell  fracture  model. 


Table  2.  Constituent  properties  for  IM7/5240-4  cubic  cell  composite 


IM7  Fiber 

5250-4  Matrix 

Composite 

E;„(GPa) 

276 

3.45 

167 

Erv,E«(GPa) 

27.6 

11.0 

V;rv,  V;,, 

0.300 

0.350 

0.320 

Vv: 

0.800 

0.510 

Gxv^  Gjtz(GPa) 

138 

1.28 

5.33 

G,,(GPa) 

7.67 

2.72 

a^(/°C) 

-0.0360x1 0’^" 

0.372x10-^ 

5.04x10-'’ 

46.8x10-'’ 

24.3x10-'’ 

To  calculate  the  Mode  1  fracture  toughness  for  the  unidirectional  composite, 
periodic  boundary  conditions  are  applied  on  the  ;c={0,  L)  and>^{0,  L]  surfaces.  The  z- 
displaeement  is  constrained  on  the  z=-L  surface,  and  the  z-displaeement  is  incrementally 
prescribed  on  the  z=L  surface.  Denote  the  applied  displacement  Uz  and  consider  a 
function  T(Az/)  defined  parametrically  with  a  parameter  w^as  follows: 

...  r  T  =  (cLzC^^z) 

■  Iau  =  u,  -  H[iC-Ha)),  +  Stnermal]  (  ^2) 

where  (o')zz  is  the  zz-direetion  volume  average  stress  in  the  unit  cell,  H  is  the  height  of 
the  unit  cell  (2L  in  the  cubic  ease),  and  Zthemai  is  the  residual  strain  due  to  curing.  The 
displacement  gap  (Au)  is  the  difference  between  the  total  displacement  and  the 
displacement  due  to  linear  elastic  mechanics,  i.e.,  the  displacement  gap  is  the 
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displacement  due  to  fracture  of  the  composite.  It  is  clear  from  equation  (22)  that  the 
composite  stiffness  and  thermal  parameters  must  be  known  a  priori  to  calculate  the 
displacement  gap.  The  unidirectional  properties  are  computed  according  to  [1],  and  thus 
the  effective  stiffness  matrix  (C^)  components  can  be  determined  from  the  auxiliary 
problems  for  the  RVE  according  to  the  formula 

efyw  =  i^/gfC/ymnCj')  +  eln)]dy,  (23) 


where  Q  represents  the  unit  cell,  and  the  w  functions  are  g-periodic  solutions  to  the  basis 
unit  strains  solving  the  six  equations 

div[C(y)(e(w^^)  +  e^^)}  =  0.  (24) 


Additionally,  the  effective  composite  thermal  expansion  coefficients  can  be  determined 
by  an  additional  computation  using  the  six  basis  solutions  namely 

Hfj  =  jCmnopCy)  +  ^op)  ■  emnO')]  dy.  (25) 

In  this  case  the  auxiliary  RVE  problem  is 


div{C(y)[e(,-n)  -  eCy)]}  =  0, 


(26) 


where  e(y)  is  the  inelastic  strain  due  to  the  mismatch  of  constituent  thermal  expansion 
coefficients.  Here,  rj  is  Q-periodic  and  e(r]')  gives  the  thermal  strains.  Figure  7  illustrates 
the  results  of  the  cubic  cell  model  fracture  simulation.  The  strength  of  the  composite  is 
predicted  by  the  maximum  of  the  volume  average  stress  value  or  maxyiugi]j{T'(Au)}.  The 
fracture  toughness  of  the  composite  is  calculated  by  integrating  the  volume  average  stress 
over  the  real  number  line  with  respect  to  the  displacement  gap  (the  area  under  the  curve 
in  Figure  7),  or 


C;c=  4nAu)d(Au)  (27) 

The  predicted  composite  strength  for  the  cubic  model  is  77  MPa,  and  the  predicted 
fracture  toughness  is  177  J/m^.  The  experimentally  determined  composite  strength  and 
fracture  toughness  values  are  66  MPa  and  225  J/m^,  respectively.  Simulations  were 
performed  at  both  the  stress-free  temperature  and  room  temperature,  with  no  difference 
being  realized  by  including  the  thermal  prestress  in  the  computation.  The  AT  for  the 
room  temperature  simulations  was  -150°C.  The  result  of  this  model  yielded  no  ability  to 
predict  realistic  strength  and  fracture  toughness  of  the  composite. 

We  concluded  a  more  realistic  microstructure  and  crack  path  is  necessary  to 
move  this  solution  toward  the  experimental  values. 
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Figure  7.  Traction-displacement  gap  plot  for  cubic  unit  cell. 


b.  Hexagonal  Unit  Cell 

After  the  cubic  unit  cell  failed  to  produce  realistic  strength  and  fracture  toughness 
predictions,  a  slightly  more  complicated  two-fiber  (total)  hexagonally  packed  cell  model 
was  employed.  The  dimensions  of  this  rectangular  prism  cell  were  lVs  x  lVs  x  L 
(jfxyxz),  with  the  origin  of  the  coordinate  system  shown  in  Figure  8.  The  location  of  the 
cohesive  zone  (fracture  surface  for  this  model)  is  again  indicated  by  the  red  arrows.  As  in 
the  cubic  model,  periodic  boundary  conditions  were  applied  on  the  x={0,  lVS}  and  y={0, 
LVS}  surfaces,  while  the  z-displacement  was  constrained  on  the  z=0  surface,  and  the  z- 
displacement  was  incrementally  prescribed  on  the  z=L  surface.  The  input  properties  for 
the  fiber,  matrix,  and  fiber/matrix  interface  domains  were  identical  to  those  used  in  the 
cubic  unit  cell  section  (see  Table  2),  however  the  ply  level  stiffness  and  thermal 
expansion  properties  varied  slightly  from  the  cubic  model  (see  Table  3).  These 
differences  are  a  result  of  the  different  fiber  microstructure  and  are  not  a  reflection  upon 
the  resolution  of  the  simulation.  The  traction-displacement  gap  simulation  results  for  the 
hexagonal  RVE  are  given  in  Figure  9.  Simulations  were  performed  at  the  stress-free 
temperature  and  at  room  temperature,  with  no  difference  in  the  prediction  realized  by 
including  the  thermal  prestress.  The  composite  strength  and  mode-I  fracture  toughness 
predicted  by  this  model  are  91.8  MPa  and  222  J/m^,  respectively.  In  this  case  the 
predicted  mode-I  fracture  toughness  is  actually  within  2%  of  the  experimental  value,  225 
J/m^.  However,  the  predicted  composite  strength  moved  in  the  wrong  direction  and  is 
now  39%  higher  than  the  experimental  value,  66  MPa. 
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Table  3.  Ply  level  properties  for  the  IM7/5250-4  hexagonal  cell  composite 


Composite 

Exx:  (GPa) 

167 

Evv9  Ej2  (GPa) 

9.70 

^xz 

0.317 

-  

0.557 

Gxv?  (GPa) 

4.97 

Gv,(GPa) 

3.12 

ax.  (/°C) 

0.369x10-'’ 

^zz 

24.8x10-'’ 

Displacement  Gap 

Figure  9.  Traction-displacement  gap  plot  for  hexagonal  unit  cell. 
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Encouraged  by  the  fracture  toughness  result  on  the  hexagonal  cell,  we  observed 
that  our  fracture  surface  was  not  particularly  realistic  for  such  a  fiber  microstructure  and 
could  thus  be  adversely  affecting  our  predicted  strength  value.  We  investigated  strength 
and  fracture  toughness  according  to  the  total  crack  length.  The  influence  of  the  intra¬ 
phase  or  matrix  crack  path  on  the  predicted  values  of  the  strength  and  fracture  toughness 
was  investigated  next.  The  fiber  packing,  material,  and  interfacial  properties  arc  identical 
to  the  previous  hexagonal  cell  model.  The  configuration  of  the  fracture  surface  is  shown 
in  Figure  10.  Here,  cohesive  zone  interfaces  are  included  on  all  fiber/matrix  interfaces. 
Figure  10. a,  and  two  mesh  independent  cracks  (MlCs  -  previously  described)  are  used  to 
connect  the  interfacial  cohesive  zones  (see  Figure  lO.b)  to  create  a  contiguous  composite 
fracture  surface. 


(a) 


Figure  10.  Configuration  of  the  fracture  surface  in  the  hexagonal  unit  eell. 

The  location  of  the  MlCs  eonnecting  the  interfacial  cohesive  zones  was  varied  to 
change  the  crack  length.  Figure  1 1  shows  a  schematic  of  the  different  MlCs  and  includes 
their  relative  crack  lengths.  Only  the  top  half  of  the  hcxagonally  packed  cell  is  included 
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in  this  schematic.  The  angles  that  are  used  to  describe  the  MICs  are  measured  from  the 
horizontal  through  the  center  point  of  the  central  fiber.  For,  example,  the  60°  MIC 
originates  at  the  central  fiber  fiber/matrix  interface  at  a  60°  angle  measured  from  the 
horizontal  of  the  central  fiber  and  terminates  at  the  upper  right  fiber  fiber/matrix  interface 
at  a  240°  measured  from  the  horizontal  of  the  center  of  the  upper  right  fiber.  A 
symmetric  MIC  was  used  to  connect  the  central  fiber  interfacial  cohesive  zone  to  the 
upper  left  fiber  interfacial  cohesive  zone.  The  65°  crack  is  essentially  the  crack  that  is 
tangent  to  both  fibers  and  thus  produces  the  shortest  overall  crack  length. 


7.000000 


Crack  Path 

Crack  Length  (um) 

Original 

13.500982 

20" 

14.907036 

25" 

13.710899 

30" 

12.766051 

35" 

12.140393 

40" 

11.783425 

50" 

11.594759 

55" 

11.499033 

60" 

11.433152 

65" 

11.426480 

Figure  1 1.  Schematic  of  mesh  independent  cracks  connecting  interfacial  cohesive  zones. 

The  traction-displacement  gap  results  of  the  MIC  fracture  simulations  for  various 
crack  path  angles  are  shown  in  Figure  12.  It  is  clear  that  both  composite  strength  and 
composite  fracture  toughness  are  affected  by  crack  length/crack  path.  Note  that  the  20° 
and  25°  crack  paths  are  longer  than  the  original  crack  path  and  have  larger  strength  and 
fracture  toughness,  while  all  crack  paths  shorter  than  the  original  predict  decreased 
strength  and  fracture  toughness.  In  general,  strength  and  fracture  toughness  appear  to  be 
linearly  related  to  the  crack  length  for  this  microgeometry  and  the  set  of  MIC  definitions. 
Minimizing  the  crack  length  decreased  the  predicted  composite  strength  from  9 1 .8  MPa 
(original  crack)  to  75.2  MPa  (65°  crack),  which  is  now  only  14%  too  large.  While  this  is 
a  massive  improvement  over  the  original  39%  error,  this  result  is  not  satisfactory.  The 
main  problem  is  that  by  minimizing  the  crack  length,  the  fracture  toughness  prediction 
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also  decreased  and  now  has  a  20%  error.  These  concerns  will  be  addressed  in  the  Future 
Directions  section  and  are  attributed  to  the  necessity  of  looking  at  more  complex  RVE. 


Volume  Average  Stress  vs  Displacement  Gap 

1.2000E-04 

l.lOOOE-04 


0.0000E4O0  5.0000E-01  l.OOOOE+00  1.5000C+00  2.0000E+00  2.5000C400  3.0000E+00  ISOOOE+OO  4.0000E+00  4.5000E+00  5.0000E+00 

Figure  12.  Traction-displacement  gap  plot  for  crack  path  experiment. 


c.  Application  to  Thermal  Oxidation 

Although  further  work  is  needed  to  improve  and  validate  the  homogenization  of 
the  fracture  toughness  and  strength  based  on  RVE  fracture  simulation,  a  strong  need  for 
such  methodology  prompted  an  early  application  of  the  proposed  technique.  This  section 
applies  the  computational  machinery  discussed  above  to  a  thermal  oxidation  problem  for 
the  AFRL/RXBC  Hybrid  Materials  for  Extreme  Environments  Program.  An  in-depth 
discussion  of  the  harsh  service  environments  of  high-temperature  composite  materials  is 
beyond  the  scope  of  this  report,  however,  references  [27-28]  describe  some  recent  micro¬ 
mechanical  work  that  has  been  performed  relative  to  these  environments  and  will  provide 
a  good  background.  The  main  issue  is  that  at  high  temperatures,  the  resin  system  will 
oxidize,  altering  the  mechanical  properties  of  the  composite.  From  a  mechanics 
standpoint  the  oxidized  resin  becomes  slightly  stiffer  and  most  importantly  significantly 
shrinks.  This  shrinkage  causes  effective  shrinkage  of  the  composite  and  creates 
significant  stresses.  These  stresses,  in  turn,  lead  to  crack  initiation,  which  provides 


■Original 
>20  Deg  Crack 
>25  Deg  Crack 
30  Deg  Crack 
35  Deg  Crack 
40  Deg  Crack 
-45  Deg  Crack 
50  Deg  Crack 
•55  Deg  Crack 
>60  Deg  Crack 
■65  Deg  Crack 
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additional  pathways  for  oxidation,  etc.  To  evaluate  the  cracking  stresses  and  simulate  the 
crack  propagation  on  the  ply  level  effective  composite  properties,  both  Group  I  and 
Group  11  properties  are  required.  In  the  case  of  traditional  composites,  as  was  discussed  in 
the  comments  pertaining  to  properties  in  Table  2,  they  can  be  measured  directly.  A 
completely  different  situation  is  encountered  in  the  problem  at  hand.  Obtaining  a 
uniformly  oxidized  state  in  a  sufficiently  large  composite  sample  without  significant 
amount  of  cracking  is  not  possible.  Only  the  neat  resin  sample  can  be  tested  in  this 
condition.  Thus  application  of  multi-scale  analysis  may  be  the  only  approach  to  solving 
coupled  thermal  oxidation  cracking  problems  in  composite  and  hybrid  materials.  In  the 
present  section  only  effective  property  predictions  will  be  addressed.  Both  Group  1  and  II 
effective  properties  for  G3-500/PMR15  composite  material  system  in  both  the  oxidized 
and  nonoxidized  regions  were  computed  as  discussed  below.  The  hexagonal  fiber 
packing  geometry  was  used,  with  a  fiber  diameter  of  7pm  and  a  fiber  volume  fraction  of 
51%.  The  65°  (tangent)  crack  path  is  used  for  the  MICs.  The  nonoxidized  and  oxidized 
properties  of  the  G3-500  fiber,  PMR15  resin,  and  resulting  composites  are  listed  in  Table 
4.  The  fiber  and  resin  properties  were  given  as  inputs,  and  all  composite  properties  were 
calculated  according  to  reference  [1]. 


Table  4.  Nonoxidized  and  oxidized  properties  of  the  G3-500  fiber,  PMR15  resin,  and 

composite 


G3- 

500 

Fiber 

Non 

Oxidized 

PMR15 

Resin 

Non 

Oxidized 

Composite 

Oxidized 

PMR15 

Resin 

Oxidized 

Composite 

E^(GPa) 

237 

3.39 

123 

4.09 

123 

E^E„(GPa) 

23.7 

7.64 

8.79 

^xz 

0.300 

0.330 

0.313 

0.330 

0.313 

0.800 

0.523 

0.536 

G^,G..(GPa) 

119 

1.27 

3.81 

1.54 

4.59 

G..(GPa) 

6.59 

2.50 

2.86 

Axial  (aiix) 
Shrinkage  (/h) 

0.00 

7.6xl0'‘’xt 

1.06xl0''xt 

2.20x10'^  + 
7.6xl0'^xt 

3.68x10'^  + 
1.27xl0'''xt 

Transverse  (yz) 
Shrinkage(/h) 

0.00 

4.24xl0-^xt 

1.23x10'^  + 
4.24xl0'^xt 

The  fracture  simulation  was  carried  out  for  both  the  nonoxidized  and  oxidized 
composite  materials  at  time  values  of  200,  600,  and  1000  hours.  The  experimentally 
determined  strength  and  Gic  data  for  the  nonoxidized/oxidized  resin  were  40/20  MPa  and 
500/250  J/m^,  respectively.  Table  5  reports  the  G3-500/PMR15  strength  and  Mode  I 
fracture  toughness  predictions  for  the  nonoxidized  and  oxidized  composite.  The  “Non” 
refers  to  the  nonoxidized  composite,  and  the  input  fiber/matrix  interfacial  strength  and 
fracture  toughness  are  equal  to  the  respective  nonoxidized  resin  values  presented  above 
(40  MPa  and  500  J/m^).  “0x1”  refers  to  the  oxidized  composite  in  which  the  fiber/matrix 
interfacial  strength  and  fracture  toughness  values  are  equal  to  the  respective  values  of  the 
oxidized  resin  (20  MPa  and  250  J/m^).  Finally,  “0x2”  refers  to  the  oxidized  composite  in 
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which  the  fiber/matrix  interfacial  strength  and  fracture  toughness  values  arc  equal  to  !4 
the  oxidized  resin  values  (5  MPa  and  62.5 

Table  5.  G3-500/PMR15  transverse  strength  and  Gic  predictions 


Non, 

200h 

Non, 

600h 

Non, 

lOOOh 

0x1, 

200h 

0x1, 

600h 

0x1, 

lOOOh 

0x2, 

200h 

0x2, 

600h 

0x2, 

lOOOh 

Strength 

(MPa) 

42.4 

42.4 

42.3 

21.2 

21.1 

20.9 

11.8 

11.7 

11.6 

GIc 

(J/m^) 

559 

556 

557 

279 

278 

277 

153 

153 

152 

The  0x2  simulation  was  used  to  test  sensitivity  to  the  fiber/matrix  intcrfacial 
strength  and  fracture  toughness  parameters.  In  the  real  composite,  the  degradation  of  the 
interfacial  properties  due  to  oxidation  is  at  this  time  unclear  at  best.  These  property 
predictions  were  passed  to  a  laminate-level  analysis  to  predict  crack  growth  and  oxidation 
layer  growth  within  the  laminate.  The  laminate-level  analysis  is  beyond  the  scope  of  this 
report.  This  multi-scale  materials  problem,  however,  represents  a  scenario  which  will 
occur  with  increasing  frequency  in  the  years  to  come  with  emerging  material  systems. 
Namely,  some  basic  information  about  the  constituent  materials  will  be  available,  but 
little  or  no  composite-level  data  will  be  available.  Computational  materials  science  will 
deliver  predictive  guidance  on  how  to  tailor  the  microstrueture  seleet  materials  to  enable 
the  revolutionary  concept  of  “the  right  material  in  the  right  spot  for  optimized 
performance.” 

d.  Future  Directions 

We  showed  that  the  microstructure  and  crack  path  play  large  roles  in  the  strength 
and  fracture  toughness  predictions.  Indeed,  wc  found  that  in  the  hexagonal  RVE,  i.c.  in  a 
uniformly  distributed  hexagonal  fiber  array,  the  crack  path  was  seemingly  not  short 
enough  to  reach  the  experimentally-determined  strength,  and  the  crack  path  was  not  long 
enough  to  maintain  the  fracture  toughness  near  the  experimental  value.  This  situation 
appears  to  be  a  dichotomy,  except  that  other  factors  in  the  microstructure  can  and  do 
affect  the  strength  and  fracture  toughness  of  the  composite.  Consider,  for  instance  a 
randomly  distributed  array  of  20  fibers,  one  realization  of  whieh  is  shown  in  Figure  13. 
Note  the  irregularity  of  the  fiber  spacing.  The  proximity  of  fibers  to  each  other  causes 
increased  stress  concentrations  in  the  matrix  domain  that  will  begin  to  fracture  the 
composite  at  much  lower  load  states  than  the  regularly-spaced  cubic  and  hexagonal 
models  presented  above.  Initiating  cracks  at  lower  load  steps  will  decrease  the  strength 
prediction  for  a  given  simulation.  As  for  the  fracture  toughness,  it  is  easy  to  see  from 
Figure  13  that  the  crack  path  will  be  more  torturous  than  the  tangent  crack  path  in  the 
hexagonal  cell,  and  thus  the  crack  length  will,  in  general,  be  longer.  Decreasing  the 
strength  while  maintaining  the  eraek  length  should  produee  predictions  that  more  eloscly 
match  the  experimental  values.  In  order  to  use  a  random  fiber  packing,  one  must  first 
verify  that  the  realization  of  the  microstructure  produces  ply-level  stiffness  properties 
consistent  with  experimental  values.  It  is  expected  that,  for  a  small  numbers  of  fibers,  the 
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Figure  13.  Sehematie  of  20  randomly  distributed  fibers. 

meehanieal  properties  will  vary  greatly  with  the  realization.  This  dependenee  should 
deerease  as  the  number  of  fibers  in  the  model  inereases.  Very  limited  work  has  been 
done  on  investigation  of  random  fiber  paeking  and  RVE  sizes  on  effective  stiffness 
properties  [6,  29],  and  the  effective  fracture  property  prediction  is  not  addressed  at  all  at 
present. 

Another  consideration  to  be  taken  into  account  is  of  a  more  fundamental  nature.  It 
concerns  the  physics  of  small  scale  cracking.  Neither  the  classical  Griffith  fracture 
mechanics  nor  the  cohesive  zone  model  used  above  are  explicitly  capturing  length  scale 
speeifie  effects,  e.g.  surface  tension,  which  may  become  important  at  the  sub-miero  and 
nano  length  scale  [30].  It  is  thought  that  the  future  development  path  of  SCSAM 
framework  should  include  such  nonelassieal  eraek  models  at  a  small  length  scale. 

Prior  to  returning  to  the  ply  level  simulations  in  the  next  section,  we  would  like  to 
return  to  the  material  property  in  Table  1  and  discuss  additional  developments  needed  to 
address  the  Group  II  properties  not  considered  above,  namely  shear  deformation  related 
strength  and  Mode  II  fracture  toughness.  The  framework  for  homogenization  of  these 
eharacteristies  is  similar  to  Eqn.  (  22);  what  is  different  is  the  eraeking  pattern  which 
occurs  in  the  RVE.  The  fracture  surface  is  expected  to  be  significantly  more  complex. 
Additional  development  of  MIC  framework  is  required  to  simulate  such  fracture 
phenomena. 
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5.  Complex  Interactive  Cracking  Network  Simulation  on  the  Ply  Level 

In  this  section  we  demonstrate  the  capability  to  model  complex  failure 
phenomena  based  on  the  MIC  and  interface  delamination  framework  on  the  ply  level.  As 
outlined  in  the  previous  section,  some  of  the  properties  required  for  the  ply  level  analysis 
cannot  be  obtained  by  means  of  homogenization  at  present.  It  is,  however,  one  of  the 
critical  practical  advantages  of  SCSAM  that  the  analysis  development  at  various  scale 
levels  ean  progress  independently  and  at  all  stages  evaluated  against  experimental  data. 

Consider  a  multilayered  composite  plate  consisting  of  N  orthotropic  layers  with  in-plane 
dimensions  L  and  W  in  the  x-  and  y-dircctions,  respectively.  Let  the  thickness  of  the  plate 
be  H  (z-direction)  as  shown  in  Figure  14.  Each  ply  represents  an  orthotropic  material 
which  is  eharacterized  by  engineering  stiffness  constants  Ey,  Gy,  Vy  and  thermal 
expansion  coefficients  tty  (ij=l,2,3).  The  direction  xi  coincides  with  the  direction  of  the 
fiber  and  the  angle  0  between  the  direction  of  the  global  coordinate  x  and  the  fiber 
direction  xi  in  a  given  ply  is  called  the  ply  orientation. 


Figure  14.  Multilayered  composite  plate  and  the  coordinate  system. 


Tensile  loading  in  the  x-dircction  will  be  applied  by  incrementing  the 
displacement  Ux  at  the  edges  x=0,  L  so  that 

u‘(0,y,z)=  uJ''(0,yz)-A‘  and  u^(L,y,z)=  uj'’ (L,y,z)+A‘  (  28) 


where  A  is  a  constant  and  /-is  the  number  of  the  loading  step.  Such  incremental 
formulation  is  required  to  properly  account  for  the  thermal  curing  stresses  prior  to 
mechanical  loading.  The  displacement  field  m/  appearing  in  equation  (2)  is  computed  by 


28 


solving  a  thermal  mechanical  expansion  problem  under  boundary  conditions,  which 
simulate  free  expansion  and  only  restrict  rigid  body  motion,  i.c. 


uj^(0.0.0)=0.  uf(W.0.0)=0  and  uj’(x,y,z)=0.  (  29) 

The  incremental  loading  boundary  conditions  (  28)  are  supplemented  with 
constraint  conditions  on  the  other  displacement  components  at  the  lateral  edges  x=0  and 
L,  so  that 


Uy(0,y,z)=  Uy^(0,yz)  and  Uj(L,y,z)=  uJ^(L,y,z),  (  30) 

which  means  that  all  displacement  components  except  axial  are  fixed  at  the  edges.  The 
analysis  begins  from  an  initially  undamaged  state  (Figure  15. a).  As  a  result  of  load 
application,  matrix  cracks  will  form  in  each  ply  of  the  laminate  as  shown  in  Figure  15.b. 
Their  location  is  arbitrary  and  not  predetermined  before  the  analysis.  At  some  value  of 
the  applied  load,  the  delaminations  between  plies  will  appear  as  a  result  of  matrix 
cracking,  as  shown  in  Figure  15.c,  and  cause  disintegration  of  the  composite  laminate. 
Note  that  some  of  the  plies  may  have  fiber  orientation  0=0  coinciding  with  the  loading 
direction,  in  which  case  these  plies  will  continue  to  carry  load  after  complete 
delamination  of  the  specimen.  The  fiber  failure  mode  will  not  be  treated  in  the  present 
manuscript  and  is  deferred  to  further  work. 


Figure  15.  Damage  progression  sequenee  in  a  laminated  composite  plate  subjected  to 
tensile  loading,  (a)  initial  stage  without  damage,  (b)  matrix  cracking  stage, 
(c)  dclamination  stage,  linking  up  matrix  cracks  in  various  plies. 
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The  initial  stress  fields  (prior  to  crack  initiation)  in  the  flat  laminates  under  axial 
tension  conditions  are  highly  uniform  in  the  x-direction,  and  therefore  the  initial  crack 
insertion  locations  have  a  tendency  to  cluster  depending  on  the  last  digit  of  computer 
number  representation.  To  obtain  more  realistic  initial  crack  insertion  patterns  and  mimic 
inevitable  statistical  variability,  quasi-random  strength  properties  were  generated  across 
the  coupon  volume.  Distributions  of  random  transverse  strength  properties  yt  (tensile),  yc 
(compressive)  and  s  (shear  strength)  were  assumed  to  follow  classical  two-parameter 
Weibull  law  defined  as 


P{x)  =  1  -  exp 


(  \ 

V 

X 

f 

1 

K 

>: 


x  =  y„yc^s 


(31) 


where  Ax  and  /3x  are  seale  and  shape  parameters  of  corresponding  strength  properties.  The 
same  value  of  12  was  assumed  for  all  strength  distributions,  based  on  transverse 
tensile  strength  scaling  in  carbon  epoxy  composites.  Using  average  strength  values  X, 
shown  in  Table  1  for  the  IM7/8552  composite,  given  in  Ref  [32],  scale  parameters  were 
calculated  as 


A^=X/r{l  +  V/3^),-,  X=Y„Ye,S 

To  ensure  mesh  independence  of  generated  quasi-random  strength  characteristics 
[31],  they  were  normalized  by  the  reference  volume  Vo  =6250mm^  (typical  of  an  8-ply 
unidirectional  composite,  which  is  used  for  transverse  strength  measurement)  and 
corresponding  local  volumes  v,.  Since  failure  criterion  was  applied  in  individual 
integration  points,  a  “local  volume”  v,-  was  associated  with  each  integration  point  as  a 

product  of  the  Gauss  weight  and  Jacobian  so  that  v,.  =  V„  ,  where  INP  is  the  total 

number  of  integration  points  and  V„  is  the  volume  of  the  n-th  ply.  Note  that  the  random 
distribution  of  initial  strength  (  31)  was  only  utilized  for  crack  origination  purposes. 
There  was  no  variability  introdueed  in  the  deeohesive  law.  This  simplified  definition 
requires  modifications  for  systematic  studies  of  stochastic  strength  distribution  effects  on 
the  apparent  strength  of  composites  and  is  used  below  only  for  method  demonstration 
purposes. 

The  methodology  presented  above  was  applied  to  perform  the  fracture  analysis  of 
the  same  quasi-isotropic  laminate  as  considered  in  [32]  and  in  the  section  above  but 
without  any  predefined  matrix  craeking  patterns.  The  analysis  methodology  correctly 
predicted  both  the  sequence  of  the  delamination,  the  failure  loads  and  the  multistep  load 
drop  behavior  before  the  final  failure  that  was  experimentally  observed  in  [32].  The  load 
carried  after  the  delamination  proeess  eorresponds  to  the  load  carried  by  the  0  plies  until 
fiber  failure.  Figure  16  displays  the  z-direction  displacement  map  during  the  stochastic 
crack  insertion  simulation  for  three  simulation  load  steps.  The  first  two  load  steps  (305 
MPa  and  406  MPa)  show  accumulation  of  the  matrix  cracking  and  delamination  onset  at 
the  free  edge  and  matrix  crack  pockets.  The  third  load  step  (495  MPa)  immediately 
precedes  the  full  delamination  of  the  -45/0  interface  and  the  associated  load  drop. 
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Figure  16.  Vertical  (z-direction)  displacement  map  at  the  different  load  steps;  (a)  305 
MPa;  (b)  406  MPa;  (c)  495  MPa. 
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Appendix:  Composite  Scarf  Repair  Optimization  for  the  JSF 

This  appendix  contains  a  final  summary  of  the  composite  scarf  repair 
optimization  that  was  chiefly  performed  during  the  previous  grant  period.  Final 
experimental  data  and  optimized  strength  predictions  were  not  available  at  the  time  of  the 
previous  final  report  and  are  included  here.  For  a  comprehensive  description  of  the  work, 
please  refer  to  reference  [33],  which  also  includes  numerical  results  concerning  patch- 
cure  warpage,  overply  thickness,  and  overply  fiber  orientation  angle  which  are  not 
presented  here. 

The  scarf  repair  is  used  to  restore  strength  to  a  composite  material  after  it  has 
been  damaged.  The  repair  process  removes  the  damaged  material  and  fills  the  void  with 
an  inverted,  truncated,  conical  patch  which  is  adhesively  bonded  to  the  remaining  original 
material  (see  Figure  Al).  The  current  standard  repair  involves  at  least  one  double 
overlay  ply  (overply)  to  restore  the  maximum  amount  of  strength.  In  Figure  Al,  the 
green  area  represents  the  original  material,  the  yellow  is  the  adhesive,  the  red  is  the  scarf 
repair  patch,  and  the  blue  is  the  overply. 


Figure  Al .  Composite  scarf  repair  models  (with  and  without  the  ovcrply). 

One  problem  with  the  current  scarf  repair  system  is  the  necessary  inclusion  of  the 
overply  to  restore  strength.  It  is  clearly  seen  in  Figure  Al  that  the  overply  adds  thickness 
to  the  original  material  and  changes  the  outer  mold  line.  This  situation  adversely  affects 
the  stealth  properties  of  the  structure,  as  well  as  increases  the  weight  and  size  of  the 
repair.  One  goal  of  this  work  is  to  improve  the  strength  of  the  repair  without  the  overply 
to  alleviate  the  negative  impacts  the  overply  has  on  the  other  critical  properties  of  the 
structure. 

It  was  previously  found  [33]  that  the  nonlinear  adhesive  modeling/failure  is  a 
critical  factor  in  composite  scarf  repair  strength  prediction.  The  nonlinear  elastic 
constitutive  equations  for  the  adhesive  material  consider  only  variation  of  the  Young’s 
modulii  and  assume  constant  Poisson’s  ratio  in  accordance  with  the  experimental 
observations.  The  shear  modulus  G  is  expressed  as 

c  =  (Al) 

where  the  elastic  modulus,  Ej,  is  a  monotonic  function  of  strain  calibrated  using  cubic 
spline  interpolation  from  experimental  tensile  testing  and  is  expressed  through  the 
dilatational  strain  invariant,  J/.  Similarly,  is  the  shear  modulus  calibrated  based  on 
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KGR-1  experimental  test  results  and  is  expressed  through  the  distortional  strain  invariant 
where 


h=  £i+  £2+  £3, _  (A2) 

J2  =  -  ^2)^  -  (£3  -  £2)^  -  (£1  -  £3)^), 

and  e,  are  the  prineiple  values  of  strain.  Equation  (Al)  is  also  used  to  predict  the  failure 
of  the  adhesive.  Critical  values  of  the  strain  invariants  (J/  and  Jj)  are  backed  out  via 
simulation  using  the  experimental  testing  on  the  neat  adhesive.  After  strain  values  result 
in  invariants  exceeding  one  of  these  critical  values,  the  shear  modulus  is  set  to  a  small 
positive  value.  Reference  repair  strength  test  values  were  reproduced  by  using  such 
adhesive  constitutive  modeling  with  accuracy. 

The  time  of  such  nonlinear  solution  was  still  unfeasible  for  optimization  studies. 
Instead,  the  chosen  goal  function  minimizes  the  elastic  L'’  norm  of  the  dcviatoric  stress 
invariant  in  the  adhesive  regime.  This  approach  determines  optimal  repair  fiber 
orientations  in  a  cost-efficient  manner.  For  strength  prediction  purposes,  a  fully 
nonlinear  analysis  is  performed  for  the  optimized  layup. 

Scarfing,  repair,  and  testing  of  the  specimens  were  performed  in  the  Materials  & 
Manufacturing  Systems  Support  Division  (RXSA).  Normalized  experimental  strength 
data  are  presented  in  Figure  A2.  The  virgin  material  strength  reference  is  for  the  ASTM 
standard  25.4  cm  x  2.54  cm  coupons.  All  other  tests  used  the  large  57  cm  x  13  cm 
coupon  size.  The  virgin  material  is  used  only  as  a  reference  data  set  to  determine 
material  properties  and  calibrate  the  critical  failure  volume  parameters  [34].  Strength 
prediction  values  closely  match  the  experimental  values  in  Figure  A2.  The  results  in 
Figure  A2  confirm  that  the  optimization  increases  the  strength  of  the  repair.  Note  that  the 
orange  bar  is  the  standard  repair  wi/h  an  overply,  and  the  red  bar  is  the  optimized  repair 
without  an  ovcrply.  The  optimized  repair  without  an  ovcrply  is  10%  stronger  than  the 
standard  repair  with  an  overply.  Additionally,  the  repair  is  20%  lighter,  10%  smaller,  and 
is  flush  with  the  original  material  surface  and  thus  docs  not  affect  the  mold  line  (stealth 
properties)  of  the  structure. 

As  a  result  of  the  success  of  this  work,  the  JSF  Repair  Group  (headed  at 
NAVAIR)  is  testing  our  optimized  repair  layups  for  use  on  the  JSF  structure.  One 
challenge  the  JSF  repair  faces  is  how  the  adhesive  behaves  in  hot  and  wet  environments. 
As  the  temperature  increases  and  the  environment  moistens,  the  adhesive  becomes  much 
more  compliant  and  thus  fails  at  smaller  loads.  Since  this  optimization  methodology 
results  in  a  repair  that  minimizes  the  adhesive  stress,  it  is  a  natural  fit  for  such  an 
environment.  The  results  of  the  JSF  testing  are  currently  not  available. 
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